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i—j We investigate the use of inexact solves for interpolatory model reduction 

00 and consider associated perturbation effects on the underlying model reduc- 

tion problem. We give bounds on system perturbations induced by inexact 
solves and relate this to termination criteria for iterative solution methods. 
We show that when a Petrov-Galerkin framework is employed for the in- 
exact solves, the associated reduced order model is an exact interpolatory 

^ model for a nearby full-order system; thus demonstrating backward stability. 

We also give evidence that for "^-optimal interpolation points, interpolatory 
model reduction is robust with respect to perturbations due to inexact solves. 
Finally, we demonstrate the effecitveness of direct use of inexact solves in op- 
timal %2 approximation. The result is an effective model reduction strategy 

_£T that is applicable in realistically large-scale settings. 

Keywords: Model reduction; system order reduction; tangential 
interpolation, iterative solves, Petrov-Galerkin 
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1. Introduction 

# > 

K/ The simulation of dynamical systems constitutes a basic framework for 

the modeling and control of many complex phenomena of interest in science 
and industry. The need for ever greater model fidelity often leads to compu- 
tational tasks that make unmanageably large demands on resources. Efficient 
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model utilization becomes a critical consideration in such large-scale problem 
settings and motivates the development of strategies for model reduction. 

We consider here linear time invariant multi-input/multi-output (MIMO) 
systems that have a state space form (in the Laplace transform domain) as 

Find v(s) such that %(s) v(s) = S(s)u(s), then y(s) = G(s) v(s). (1) 

Here, u(s) and y (s) denote Laplace-transformed system inputs and outputs, 
respectively; v(s) represents the internal system state. We assume that 
G(s) G C pxn and B(s) G c nxm are analytic in the right half plane; and 
that DC(s) G c nxn is analytic and full rank throughout the right half plane. 
Solving for y(s) in terms of u(s), we obtain 

y(s) = e(s)X(s)- 1 S{s)u{s) = 3€(s)u(s). (2) 

This representation of the transfer function, 

K(s) = e{s)DC{s)- 1 'B{s), (3) 

we refer to as a generalized coprime realization. Standard first-order descrip- 
tor system realizations, with J€(s) = C (sE — A)" B for constant matrices 
E, A G R nxn , B G R nxm , and C G R pxn evidently fit this pattern with 
C(s) = C, S(s) = B, and DC(s) = sE — A. However, many dynamical sys- 
tems can be described more naturally with generalized coprime realizations. 
For example, a system that includes internal system delays as well as trans- 
mission/propagation delays in its input and output could be described with 
a model 

Ei(t) = A x(t) + Ai x({ - r S!/s ) + B u(t - r mp ), y(£) = Cx(£ - r out ) (4) 

for r sys , r mp , r out > 0, and E, A , Ai G R nxn , B G R nxm and C G R pxn . 
Taking the Laplace transformation of Q yields the transfer function 

M{s) = e(s)3C(s)^ 1 S(s) = ( e - STout C) (sE- A - e" 87 "^ Ai)^ 1 (e~ ST ^ B) , 

which has the form of (J3J). The form of (J3J) can accomodate greater gen- 
erality than this, of course, including memory convolution involving higher 
derivatives, second and higher-order polynomial differential equations, sys- 
tems described via integro-differential equations, and systems where state 
variables may be coupled through infinite dimensional subsystems (possibly 



Table 1: Examples of Generalized Coprime System Realizations 



Descriptor Systems 


C(sE — A) X B (E possibly singular) 


Delay Systems 


(e- sr »"C)(sI - A - e- aT ""Ai)- 1 (e- STmi 'B) 


Second Order Systems 


(sCi + C )(s 2 M + sG + K) ^B 


Weighted Systems 


W ( S )C(sI-A)- 1 BW i ( S ) 



modeling internal propagation or diffusion). See Table [I] for other examples 
and [6] for further discussion. 

In many applications, the state space dimension, n, is too large for efficient 
system simulation and control computation, so the cases of interest for us here 
have state space dimension vastly larger than input and output dimensions: 
n 3> m,p. See [19] for a recent collection of such benchmark problems. 

The goal is to produce a reduced system that will have approximately the 
same response (output) as the original system for any given input u(t). For 
a given reduced-order r <C n, we construct reduced order models through a 
Petrov-Galerkin approximation of (flh: Select full rank matrices V r G R nxr 
and W r G M. nxr . For any input, u(t), the reduced system output, y r (t), is 
then defined (in the Laplace transform domain) as: 

Find v(s) G Ran(V r ) such that W^ (%(s) v(s) - S(s)u(s)) = (5) 
theny r (s) =e(s)v(s) (6) 

which defines the reduced transfer function as, 

K r (s) = e r (s)JC r (s)- l 'B r (s), (7) 

where 



^rxm 



% r ( s ) = W' r %(s)V r G C rxr , S r (s) = W; 3(a) G 

and e r (s) = C(s)V r G C pxr . (8) 

2. Interpolatory Model Reduction 

Interpolatory reduced order models are designed to exactly reproduce 
certain system response components that result from inputs having specified 
frequency content and growth. The approach has been described for standard 



first-order system realizations in [T31 121 EH E] and extended to generalized 
coprime realizations in [6] . We summarize the basic elements of this approach 
below. 



A set of points {Hi\ r i=l C C and (nontrivial) direction vectors {cj} 



i=l 



C 



constitute left tangential interpolation data for the reduced model, 9£ r (s), if 
cfJ€(/ij) = cf IK r (/ij) for each i — 1, . . . , r. (9) 



Likewise, {(Jj} r j=v and associated directions {bj}^ =1 C C m , constitute right 
tangential interpolation data for the reduced model, 3€ r (s), if 



3-C(oj )bj = JKj.(cTj)bj for each j — 1, 



r. 



(10) 



Given left and right tangential interpolating data, interpolatory model re- 
duction may be implemented by first solving the linear systems: 



Find Wj such that w^ OC(fii) = c { C(/i 



find Vj such that 3C(<jj)vj 



Sfo)bi 



for i = 1, 

for j -- 



r, and 



r. 



(11) 
(12) 



We assume that the two point sets {^i}\ =1 and {(Tj} r - =1 each consist of r 
distinct points and that the vectors {vi, ••• , v r } and {wi, ••• , w r } are 
linearly independent sets. These vectors constitute "primitive bases" for the 



subspaces V r = spanjvi, 
associated matrices: 



v, r } and W r = spanjwi 



w r }. Define the 



V r 



[v 



i- 



W 



r 



Wi 



w; 



[0C(a 1 )- 1 'B(a l )b 1: 



cjG(fi r )0C(fi r 



%{a r )- l S{a r )b r ] , (13) 
(14) 



The reduced model, ^K r (s), as defined in ^ and (J8| using V r and W r from 
(13) and (14), interpolates J€(s) at the 2r points {(j>i} r i=1 and {<Jj} r j=1 , in 



respective output directions {cj} i=1 and input directions {hj} r j=1 ; that is, 
conditions ^ and (10) are satisfied. If /i^ = a k for some k then first order 
bitangential moments match as well: 

cl J<'(n k ) b k = cl ( K' r {fx k ) b fc 

Interpolation of higher order derivatives of IK(s) can be accomplished with 
similar constructions as well; see [HI [3] and references therein. 



For large-scale settings with millions of degrees of freedom, interpolatory 
model reduction has become the method of choice since it does not require 
dense matrix operations; the major computational cost lies in solving the 



(often sparse) linear systems in (11 ) and (12). This contrasts with Gramian- 

opti- 



based model reduction approaches such as balanced truncation 
mal Hankel norm approximation [T2] and singular perturbation approxima- 
tion [21] where large-scale Lyapunov equations need to be solved. Moreover, 
these computational advantages have been enhanced for standard first or- 
der state-space realizations by strategies for optimal selection of tangential 
interpolation data, see jTB]. 

2.1. Inexact Interpolatory Model Reduction 

The basic framework for interpolatory model reduction presumes that the 



key equations (11 ) and (12) may be solved exactly or nearly so, at least to an 



accuracy associated with machine precision. Direct solution methods, em- 
ploying sparse factorization strategies, for example, are capable of handling 
systems of significantly large order. However since the need for ever greater 
modeling detail and fidelity can drive system order to the order of millions, 



the use of direct solvers for the linear systems (11) and (12) often becomes 



infeasible and iterative methods must be employed that terminate with pos- 
sibly coarse approximate solutions to the linear systems. We consider and 
evaluate issues related to these approaches here. 

Suppose {vi, • • • , v r } and {wi, • • • , w r } are linearly independent sets 
in C n and define 



V r = [Vi, 



W, J 



w 



r 



w 



T 



(15) 



Wj and Vj will be viewed as approximate solutions to the linear systems (11 ) 



and (12) and accordingly we will refer to them as "inexact" solutions to (11) 



and (12). Nonetheless, unless otherwise stated, these vectors can be any 



arbitrarily chosen linearly independent vectors in C n . 

Define residuals, £ { and ry ■, corresponding to Wj and Vj, as 



%{a j )^ j 



S(^)b, 



& = X(jii) ^ - e(fii) c % and rjj 
The deviations from the corresponding exact solutions are then 



5wi 



Wi 



w, : 



X(fj,i) T £j and 5vj = Vj — Vj 



%{*,) 



Vi 



(16) 
(17) 



The resulting (inexact) basis matrices destined for use in a reduced order 
model are 

W r = W r + [Sw h ■■■ , 5w r ] (18) 

V r = V r + [«Jvi, ••■ , 6v r ]- (19) 

Define reduced order maps associated with these inexact bases: 
3C r (s) = W r T 3C(s)V r , S r (s) = WjS(s), and G r (s) = C(s)V r , (20) 
together with the associated inexact reduced order transfer function 

$C r (s) = e r (5)3C r (5)~ 1 S P (s). 

Notice that we are free to make any choice for bases for the subspaces, V r and 
W r , in defining J€ r (s); no change in the definition of ( |20 ) is necessary. As a 
practical matter, it is generally prudent to choose well conditioned bases in 
computation. 

3. Forward Error 

3.1. Interpolation Error 



Inexactness in the solution of the key linear systems (11) and (12) pro- 
duces a computed reduced order transfer function, ^K r (s) that no longer in- 
terpolates JC(s); typically, the reduced order system response will no longer 
match any component of the full order system response at any of the complex 
frequencies {//j}[ =1 and {ctj}[ =1 that have been specified. How much response 
error has been introduced at these points ? 

The particular realization taken for a transfer function can create innate 
sensitivities to perturbations associated with that representation. Define 
perturbed transfer functions, 

9C«s(s) = e{s)3C(s)-\'B(s)+S'B) and 3<se(s) = (e(s)+5e)DC{s)- 1 'B(s). 

In discussing perurbations in system response caused by 8"B and 6G at s = a, 
it is natural to introduce the following quantities: 



cond:s(3-C(o-)) 
conde (%(cr)) 



le^DC^)- 1 !!!!©^)!! 



|W( 



<j 



(a)\\\\JC(ci)- l 'B(a) 



\\n 



<J 



to be condition numbers of the transfer function response, by way of anal- 
ogy to the condition number of algebraic linear systems. (Unless otherwise 
noted, norms will always refer to the Euclidean 2-norm for vectors or the nat- 
urally induced spectral norm for matrices). It is straightforward to show that 
these quantities measure the relative sensitivity of the system with respect 
to perturbations in "B and C, respectively: 



IIKWII - v v "ll»M 

l%eW-MWIU conde(K(<7 || M!| 



ll«WII - " v "liewn 

For values of s such that 3C r (s) and JC r (s) are nonsingular, define the 
matrix-valued functions, 

9 r (s) = 3C(s)V r 3C r (s)- 1 Wj, Q r (s) = VrDCris)- 1 WjOC(s), 
9 r (s) = JC(s)\ r 6c r (s)~ 1 Wj, and Q r (s) =\ r 0C r (s)- 1 WjJC(s) (21) 

Where defined, T r (s), Q r (s), T r (s), and Q r (s) are differentiable (indeed, 
analytic) with respect to s, having derivatives that satisfy: 

$' r (s) = (l - 9 r ) %'(s)X(s)- 1 T r 
and V J (22) 

Q r (s) = Q r 0C(s)- 1 DC / (s) (l - Q r ) 

with equivalent expressions for ^P r (s) and Q' r (s). We will make a series of 
observations about properties of 3\-(s) and Q r (s) which will have immediately 
apparent parallels to properties for T r (s) and Q r (s). 

Observe first that T r = T r and Q r = Q r so both T r (s) and Q r (s) are 
skew projectors. These projectors are of interest because the pointwise error 
in the transfer function can be expressed as 

W(s) - $C r (s) =Q(s) (jCis)- 1 - V r DC r (s) _1 W^ B(s) 
=e(s)3C(s)- 1 (l - P r (s)) B(a). 

Similarly, 

M(s) - ttr(s) = 6(a) (i - Q r (s)) 3C(s)- 1 S(s) 



and 

W(s) - St r (s) = e(s) (i - Q r (s)) JC(s)- 1 ( I - iP P (s) ) S(s). 



The derivative of this last expression can be computed with the aid of ( 22 ) 
and observing X(s)- 1 9 r (s) = Q r (s)3C(s)- 1 : 

W(s) - #£(s) ~ [e(s)aC(s)- 1 ] (i - $ r (s)) S(s) (23) 

+ e( s )(i-Q r ( s ))^[3<;( s )- 1 s( s )] 

- e( S ) (i - Q P («)) ^ [^(s)- 1 ] (i - $ r ( s )) B( s ). 

We introduce the following (s-dependent) subspaces: 

<P P (s) = Ran 9 r (s) = Ran 9C(s)V r , Q r (s) = Ker (W^K(s)) 1 , 

$ r (s) = Ran y r (s) = Ran 9C(s)V r , £> r (s) = Ker (wj3C(s)) , 

<8 TO (s) = Ran 3C(s)- 1 B(s), £ p (s) = Ker (e(s)OC(s)- 1 ) 1 . 

CP r (s) maps vectors in C n onto *P r (s) along W,f and Q r maps vectors in C n 
onto V r along £l r .(s)- L . 

Given two subspaces of C n , say M. and TV, we express the proximity of one 
to the other in terms of the angle between the subspaces, Q(A4,Af) G [0, |] 
defined as 

sup inf ,, ,, = sinGfA^AO- 

xex y 6Ar ||x|| 

Q(A4,Af) is the largest canonical angle between M. and a "closest" subspace 
A/" of M having dimension equal to dim A4. Notice that if dim M < dim A4 
then <d(M,Af) = f and Q(M,Af) = if and only if M C Af. Q(M,Af) 
is asymmetrically defined with respect to KA and A/", however if dim TV = 
dimA4 then @(M,Af) = Q(Af,M). If Hm and IIv denote orthogonal 
projectors onto M. and Af, respectively, then sin Q(A4,Af) = ||(I— IIm)II/v||. 
The spectral norm of a skew projector can be expressed in terms of the 
angle between its range and cokernel [27]. In particular, 



|$ r 00|| = 111-^0011 = —^— ^- (24) 

cos9(<p r (s), W r ) 

|Q P ( S )|| = ||I - Q P ( S )|| = —J— (25) 

cos6(n r .(s), V r ) 



Theorem 3.1. Given the full-order model !K(s) = C(s)3C(s) _1 B(s) ; inter- 
polation points {o~j} C C, {fa} C C and corresponding tangential directions, 
{bj} C C m and {c^} C C p , let the inexact interpolatory reduced model 

Ji r (s) = C r (s)9C r (s) _1 S r (s) be constructed as defined in (15)-(20). The 
(tangential) interpolation error at fa and o~j is 



\!}{. r ((Tj)hj — 5C(c,)b, II , ,_„. 

11 ^- <cond' B {'K{a j )h j 



sin 8 (<£p(<Tj), W, 



\V_i\ 



\\K(<Tj)bj\\ 

|cfJ€ r (/xi) - cf JC(/xj 
llcfMCw)!! 



cosefeo.,), w r ) ll s (°j) b il 



<cond e (cf 3{(fa)) 



sinG MB m (/Xj), V r 



cose(5 r (/xi), Pr) ll c i e (/^ 



(26) 
• (27) 



1/ /ij = o-j i/ien, 
|cf JCr(/ii)bi - cJJC(fa)bi\ < 



i^n wviWwti 



maxl cosG(*p r (/ij), W r ) , cos0 (0. r (fa), V., 



(28) 



and 



cj W(jn)bi - cj M r (Mi)bi| < M 



+ 



Wi 



-=- + 7^ ~ 

\Vi\\ ll&ll 



cose(q3 r (/ii), Wr) cose(Q r (/ii), v r 



(29) 



with M = mail 



[cfeac- 1 ] 



c/.s 



[JC'SbJ 



[X-'] 



Proof: From (17), v, = OC(o-j) 1 (S(a j )b j +rf j ), which implies then that 



3C(o-j)vi = ©(o-jObj + r/j- G qj r (o-j) and (I - iPrfo)) ( s (°j) b j + fj) = °> 
which may be rearranged to obtain 



I-? r fo) Sfo) b i 



I-DVfo) )r/,. 



(30) 



Let II be the orthogonal projector taking C n onto W r = Ker i'J > r (s) ) . One 



may directly verify that I — 3V(s) = [I — II ) ( I — T r (s) ) , and 

$Lr{?j)hi - Wfajb, = -e(a J )3C(^)- 1 (i - Tr^)) Bfojb, 

= e(a J )OC(a J )- l (l-$ r (o J ))r 1j (31) 

= eMXfa)- 1 (i - fi) (i - $ r (a s j) rjy 

Now suppose r is an orthogonal projector onto £p((Tj). We have then that 
Ran(I - T) = Ker(e(a i )0<:(a i )- 1 , so that &(a j )'K(a-)- 1 ) = e(a j )DC(a i )- 1 r 
and 

fcriff^bj - 'K{a j )h j = e(a j )DC(a i )- 1 r (i - ft) (i - 9 r {aj)\ rjj. 



Taking norms, we obtain an estimate yielding (26): 



l^iria^hj - XiarfbjW <\\(l-U\T {e{a j )0C{a j )- l ) T \\ ■ ||I - iP r (o>-) II ' \\Vj 

sine (tpto-j), W r 

< iie^ox^-)- 1 !! ^ ^J- ■ ||»7 7 . 

cose(^ r (aj),Wr" 



(27) is shown similarly, noting first that 

cfe(fM) (i - Orfcii)) = -& (i - &.(/**)) ■ (32) 



Defining the orthogonal projector, H, that takes C n onto V r = Ran ( Q r (s) ) , 
one observes next I — Q r (s) — ( I — Q r (s) ) ( I — E] so that 

Hcf^^) - cfiHto)!! = llcfe(^) (i - Q r (/i,)) acG^sfrOII 

< ||£f (i - Q r (^)) (i - §) XOk)- 1 *^ 

< U % \\ ■ ||I-Q r (^)|| • || (i - s) Xbur^ini) 

sin9(5S m (/Xi), V r 

< IIDC^) -1 ©^)" 



COs6(n r (/ii), V T 



10 



When fa = di, we have 

= if (i - o r (^)) xOh)- 1 (i - $r(mj) rti 

^( /Ul )- 1 (l-y r ( / u i ))r ?i , or 

^(i-o,^))^^)- 1 ^, 

leading then to two estimates: 

\cJ'K{^)h l -c^ r ^ i )h i \ < ||e,|| • 11*7,11 • H3CO0" 1 !! ■ ||I-^(^)|| 
and 

IcfM^hi-c^Mp^bil < ||t<|| ■ II^H • UXCAii)- 1 !! ■ ||I-OrO*i)l|. 



These can be combined to yield ( 28 ) . 



The last inequality comes from using (23) with s = jjf. 



cfWGuObi - cf£t(#)b< = j- s [cjGJC- 1 ] | w (i - P r (^)) S(^)bi 

+cf e(/i.) (i - Q r .(^)) ^ pe^sbj | w 

-cf e^) (i - Q r (/i,)) ^ [x- 1 ] | w (i - £.(/*)) S^b*. 



Then from ( [30] ), (J32J), and the Cauchy-Schwarz inequality 



c^JC0* i )b i -c^JC r (Mi)b < 



< 



rf.s- 



[cfCDC- 1 ]! (l-PMjrj 



+ 



d 



^(i-Q^M-ix-^m 



+ 



$ (I - Q r ^)J - [X- 1 ] ^(1-^(^)17, 



< 



,T/o <v— 11 



T f(fex , 

ds L * J I W 



+ 



l»7il 



cos0(jQ r (/Xi), V, 



COsG(qj r (Mi), W r ) 

' f ix-^biW 

L *J l/ii 



d,s 



+ 



— in 



ds L J '« 



l»7il 



cos9(q3 r (/ii), VW r ) cose (5, (Mi). v r 



11 



which yields the conclusion. □ 



Consider the effect of solving (11) and (12) approximately with succes- 



sively increasing levels of accuracy that force the residual norms to zero, 
— > and ||£j|| —?• 0. The multiplicative behavior of the error bound (28) 
r/|| and ||£J| contrasts with the additive behavior seen in 



with respect to 



(26) and (27) and suggests some potential benefit in using the same inter- 
polation points for both left and right interpolation, i.e., choosing \ii = 0[ 
for i = 1, . . . , r. Note that this choice also forces convergent (bitangen- 
tial) derivative interpolation as shown in (29). Indeed, choosing /ij = <ji for 
i = 1, . . . , r is a necessary condition for forming "H 2 -optimal interpolatory 
reduced order models for first-order descriptor realizations, as we discuss in 
§p| (see also [IE]). Beyond this, there can be notable computational advan- 
tages in choosing \ii = 0-j, since the linear systems to be solved in (11) and 



(12) then have the same coefficient matrix; allowing one potentially to reuse 



factorizations and preconditioners. 

Certain applications require the retention of structural properties such as 
symmetry in passing from 3C to 0C r and one is compelled to choose W r = V r 
("one-sided" model reduction), so the vectors {wi, ••• , w r } might not be 
approximate solutions to (|Tlj) in the usual sense. Nonetheless, the behavior 



of the interpolation error is still governed by (26) and (27). We explore this 
in the following numerical example. 



We illustrate the character of the results given in Theorem 3.1 , bounding 
the response error at the nominal interpolation points caused by inexact 



solves in (11 ) and (12). To this end, we consider a delay differential equation 
of the form introduced in (El) taking n = 2000, m — p — 1 and r mp = r out = 0. 
The coefficient matrices for the full order model in Q were taken from [BJ. 
We construct multiple reduced models all of order r = 3 , solving (11) and 



(12) with different levels of accuracy. We chose three logarithmically spaced 
values, G\ = 0.001, <Ji = 0.0316, cr^ = 1.0, and fixed them as interpolation 



points. We then obtained approximate solutions of varying accuracy to (11) 



and (12) in a manner described in more detail below, assembled the inexact 
interpolation basis matrices, V r and W r , and obtained reduced models of 
order r = 3 having the same internal delay structure as the original system: 

$Cr(s) = e r .(s)3C r .(s)- 1 S r (s) 

= CV r (sW^EVr - W^A V r - e- ST °y° W?AiV r V W^B 



We considered both the usual "two-sided" model reduction process that in- 



12 



volves approximate solution of both ( 11 ) and p2| and the "one-sided" process 



that involves approximate solutions only to (121) to generate V r and then as- 
signing W r = V r . Linear systems were solved with GMRES terminating 
with a final relative residual below a uniform tolerance denoted by e. 

We generated reduced order models in this way, varying the relative resid- 
ual tolerance e from 10 _1 down to 10~ 8 . Figure [l] below shows the resulting 
interpolation errors |!K(o~i) — !K r (o"i)| and bounds from equations (26) and 



(28) for one-sided and two-sided cases, respectively, as e varies. Observe that 



the bounds in Theorem 3^ predict the convergence behavior of the true er- 
ror quite well; the rates (slopes) are matched almost exactly. Note also that 
the interpolation error decays much faster for two-sided reduction than for 
one-sided reduction Indeed, the ratio of the two errors is close to e, i.e., for a 
given tolerance e, the interpolation error for two-sided reduction is approxi- 
mately e times smaller than the interpolation error for one-sided reduction. 



10""' 



Interpolation error vs the upper bound at a 




O Error- 1 sided 

■ ©- Bound-1sided 
-x — Error-2sided 

■ x - Bound-2sided 



10 



10"' 



10 



10 



10 10 10 10"' 



Figure 1: Behavior of interpolation error and upper bounds vs e 

Analogous results regarding behavior of the bounds and interpolation 
error are observed at Oi and a^ and so are omitted for brevity. 

3.2. Global Error Bounds 

Thus far we have focussed on the extent to which interpolation properties 
are lost in the computed reduced models when inexact solves are introduced 
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into the process, considering in effect local error bounds. Clearly, it is impor- 
tant to understand the effect of inexact solves on the overall global quality of 
the reduced order model. There are two commonly used measures for close- 
ness of two conforming dynamical systems (i.e., those with the same input 
and output dimensions): 

1 f°° 
the %2-norm: ||3~C _ S||w 2 = — / \\^K(tco) — S(tco)\\ 2 F dco 

2tt J_ oq 

the "Hoo-norm: ||!K — S||"Hoo = max || 9C(zu;) — S(«w)|| 2 . 

Since reduced models are completely determined by the subspaces, V r and 
W r , as shown in rt8|), we first evaluate (in Theorem 3.2) how much inexact 



interpolatory subspaces, V r and W r , can deviate from the corresponding 
true subspaces, V r and W r , as a result of inexact solves. The effect of this 
deviation on the resulting model reduction (forward) error will be shown in 



Theorem 3J3 In this way, we are able to connect model reduction error 
to observable quantities that are associated with inexact solves, such as the 
relative stopping criterion e. 

Theorem 3.2. Let the columns o/V r and\ r be exact and approximate solu- 
tions to (12) and the columns ofW r andW r be exact and approximate solu- 
tions to (11). Suppose approximate solutions are computed to a relative resid- 
ual tolerance of e > 0, so that ||?7 i || < £:||23(cr,;)bj|| and ||£J < e \\C(jj,i) T Ci\\ , 
where the residuals rj i and £ t are defined in (!(% )■ 

Denoting the associated subspaces as V r , V r , W r and W r then 



sine(V r , V r )< — ^- (33) 

?min(V r D„) 



sine(W r , W r ) < e -£- (34) 

where Y) v and T) w are diagonal scaling matrices defined as 

D, = diag{(\\OC(a 1 )- 1 \\\\'B(a l )b 1 \\)-\ . . . , {{{Xia^l HS^M)" 1 ) and 

n w = ^^((iiacCMir 1 !! lie^ifciH)- 1 , ..., (wn^r'W \\e^ r f e r \\r l ) 

and ^min(M) denotes the smallest singular value of the matrix M. 
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PROOF: We prove (33). The proof of (34) is similar. 



Write V r = V r + E with E = [K^i)" 1 ^, . . . , K(a r )- 1 r 7r ]. Then 

-^ ||v — v|| 

sm0(V r , V r ) = max min 



max mm 



vev r vev r ||v|| 

Ya=i z i OC(a i y 1 'B(a i )hi - Ya=i x fii 



X i Z i || 2^/1 = 1 x i' v i\\ 

ELi(^ - x i )X(a i )- 1 'B(ai)b i - xtfC^)' 1 ^ 



: max mm 



i=l X i V i 



ll£Li*i 3C W~ 1, fcll l|Ex|| ||EDx|| 
< max r ^rj = max — ^ = max — ^ 

x * \\2^i=l X i V i\\ x ||V r x|| x ||V r Dx|| 

where D = diag(<ii, . . . , d r ) is a diagonal matrix with positive diagonal 
entries, di > 0, that are fixed but for the moment unspecified. 
Note that 



||EDx|| <||ED||||x|| < v/r||x|| max((i i ||3C(ai) _1 r7 l 

i 

< Vr||x|| m a x(d l \\%(a l )- 1 \\\\rj l \\) 

Thus we have, 

max; (dj||DC(o-j) _1 || ||ry) __ m.ayL i (d i \\'K(<T i )- l \\\\'q i 



sinG(V r , V r ) < Vr 



min x (j|V r Dx||/||x||) * ?min(V r D) 

(35) 
This bound is valid for any choice of diagonal scalings, D, so we can min- 



imize the right hand side of (35) with respect to d\, . . . , d r . The Column 
Equilibration Theorem of van der Sluis [28] asserts that the optimal choice of 
di, . . . , d r is such that d;||!K((Xj) _1 || \\rii\\ = C, independent of i — 1, . . . , r. 
If inexact solves terminate with residuals satisfying HryJI ~ e ||25(crj)bj|| then 
we may take C = e and di = (||3C(o"j) _1 || ||S(crj)bj||) to achieve the best 



bound possible with the information given. This leads to (33).D 



As a practical matter, the column scalings used in (33) and (34) will 
not be computationally feasible in realistic settings. If instead we scale 
the columns of V r and W r to have unit norm (cheap !) - taking D^ = 

diag(l/||vi||, ..., l/||vy||) andD„ = diag(l/||wi||. . , l/||w r ||), the bound 
.es to 

sin6(V r , V r ) < maxK 2 (^(c^i), Vj) 
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for (33) degrades to 



where k 2 (3C(ctj), Vj 



ix^)- 1 !!!!®^)^! 



> 1 is the condition number 



of the linear system (12). A similar expression holds for sin0(W r , W r ). In 



many cases, these condition numbers have only modest magnitude and the 



bounds (33) and (34) remain descriptive. 



Theorem 3.3. Let the columns of V r and V r be exact and approximate 
solutions to 1 12) and the columns ofW r andW r be exact and approximate 
solutions to (11). Let the associated subspaces be denoted as V r , V r , W r and 



W r and the associated reduced order systems be denoted as !K r (s) (exact) and 
IK r (s) (inexact). Then 



\n_r n/ 



r\\Hc 






\H C 



+ 3"C r 



\U C 



< M max (sine(V r , V r ), sinG(>V r , W r 



where 



M = 2 max 



max^gjgi conde(9£ r (iuj)) max^gjgi condsCK r (iuj)) 



min^^ cos e(£l r (ico), V r ) ' min o;eM cose (^M' W r 



and 



cond&(fK. r {s)) 
conde(J€ r (s)) 



lejsWJs^wiwwBis 



\G(s) 



\\H r (s)\\ 

WJKr(s)- 1 3Js) 



\K r (s) 



PROOF: Note that for all s G C for which lK r and !K r are both analytic, 
||M r («) - Mr (a) || = ||6(s) (VrXris^W^ - V r 3C r (s) _1 W^ 3(s)|| 

= ||e(s) (Qr(s) - Q r (s)) x^y 1 ®^ 

= ||e(s) (Yl - Q r (s)) Q r (s) - Q r (s) (I - Q r (s))) OCis)- 1 ^) 
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So, 
||W r (s) 



HJs)\\ < \\e(s) 



< \\e(s) i 



I-Q r (s)jO r (s)3C(s)- 1 S(s)|| 
+ ||e(s)Q r (s) (I - Q r (s)) K^)" 1 ©^) 
O r (s))Q r (s)3C(s)- 1 S(s)|| 



+ ||e( s )3C( s )- 1 y r ( s )(i-y r (s))S( s )|| 
< ||e(s) (i - o r (*)) (i - s) sq,.^)^)- 1 ©^)!! 

+ ||e(s)3C(s)- 1 5' r (s) n (i - n) (i - a> r (a)) S( 



< 116(a) 



< 116(a) 



I -Or M 






|Q r (s)3e(s) _1 B(a 



+ ||e(s)3C( s )- 1 ^^(s)|| ||n (i - n) || ||i - av(s)|| ||®( a ) 



sine(V r , V r ) 

cose(d r (s), v r ) 



|Q r (s)DC(s) _1 S(s) 



+ ||e(s)DC(s)- 1 y r (s) 



sinO(W r , VW r ) 
cos0(*p r (s), W r 



IBOO 



< cond e (J€ r (s)) ^|ff^ ||:K,.(.s) 



+ pC.(s 



cos6(n r (s), V r 

sin9(VV r , W r ) 



cos0(^ r (s), W r 
Maximizing over s = 100 with w6R gives 

sin0(V r , V r ) 



cond-B CK r (s)) 



9i r — IKrll'Hoo — maxcondg(!K r (?a;)) — 

ojgR min^gjj cos6(£2 r (s), V r ) 






He 



, ,£., sin9(W r , W r ) 

+ max cond-s -H r U<^ ) ) ^ .^, , . — — — - 

^eR V V ; Vin weM cose(^ r (s), W r ) 






rto 



which leads immediately to the conclusion. □ 



3.3. Illustrative examples 

The process to be modeled arises in cooling within a rolling mill and is 
modeled as boundary control of a two dimensional heat equation. A finite 
element discretization results in a descriptor system of the form 

Ex(t) = Ax(t) + Bu(t), y(t) = Cx(t). 
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where A, E e r5177x5177 ; b e m 5177x7 7 C G M 6x5177 . For simplicity, we focus 
on a SISO full-order subsystem that relates the sixth input to the second 
output. For details regarding the modeling, discretization, optimal control 
design, and model reduction, see 0E]. 

We show the results of interpolatory model reduction using an ad hoc 
choice of interpolation points: 6 logarithmically spaced points between 10°' 5 
and 10; and an %2-optimal choice of interpolation points obtained by the 
method of [16]. For each case, we reduce the system order to r = 6 using 
first exact interpolatory model reduction (i.e., the linear systems are solved 
directly) and then with inexact model reduction with varying choices of ter- 
mination criteria. The resulting reduced-order models are denoted by !K r (s) 
and !K r (s), respectively. To see the effect of the choice of interpolation points 
on the underlying model reduction problem, we vary the relative residual ter- 
mination tolerance, e between 10 -1 and 10~ 10 and show how quickly !K r (s) 
converges to 3-C r (s) for both the ad hoc selection and the % 2 -optimal selec- 
tion of interpolation points. Table [2] shows the relative H^ error between 
J€ r (s) and IK r (s) as e decreases. For the T^-optimal choice of interpolation 
points, < K. r {s) converges to < K. r {s) as e decreases, for the ad hoc choice of 
points, there is almost no improvement in accuracy until e = 1 x 10~ 6 . 



e 


"^-optimal {<jj} 


ad hoc {o~i} 


KT 1 


7.22 x 10" 1 


5.05 x 10- 1 


io- 2 


2.00 x IO" 1 


1.64 x IO" 1 


IO" 3 


4.27 x IO" 2 


4.11 x IO" 1 


10" 4 


1.07 x IO" 2 


2.38 x IO" 1 


10~ 5 


2.76 x IO" 4 


5.62 x IO" 1 


IO" 6 


2.56 x IO" 5 


2.13 x IO" 2 


io- 7 


2.91 x IO" 6 


3.52 x 10~ 3 


io- 8 


1.51 x IO" 7 


6.18 x IO" 5 


io- 9 


2.07 x IO" 8 


1.76 x IO" 5 


io- 10 


2.17 x IO" 9 


5.15 x IO" 6 



Table 2: The relative error 



•/If J\,p 



\jv r 



as e varies 



\H, 



The behavior exhibited in Table [2] becomes clearer once we inspect the 
subspace angles between the exact interpolatory subspaces V r , W r and the in- 



exact ones V r and W r . Table |3j shows the sine of the angle between the exact 
and inexact interpolatory subspaces as e varies. While the gap decreases sig- 
nificantly as e decreases for an %2-optimal selection of interpolation points, 
there is a much smaller improvement in the gap with respect to e for an ad 
hoc choice of points. This behavior will be re-visited in more detail in §4.2 
revealing that the "^-optimal (or good) interpolation points are expected to 
produce reduced order models that are more robust with respect to pertur- 
bations due to inexact solves. 





sinG(V 


r,V r ) 


sin6(>V 


r,W r ) 


e 


%2-optimal {a.j} 


ad hoc {a.j} 


%2-optimal {a,i\ 


ad hoc {(Ji} 


10- 1 


9.85 x 1CT 1 


9.99 x IO" 1 


9.99 x KT 1 


9.99 x 10" 1 


lO" 2 


1.99 x KT 1 


9.99 x IO" 1 


9.97 x 10" 1 


9.93 x 10~ 4 


1(T 3 


2.36 x IO" 2 


9.99 x 10" 1 


4.87 x KT 1 


9.83 x IO" 1 


10" 4 


4.39 x 1(T 3 


9.60 x 10" 1 


6.38 x 1(T 2 


9.99 x 10" 1 


10" 5 


2.72 x 1(T 4 


5.80 x IO" 1 


7.09 x IO" 3 


7.20 x IO" 1 


1(T 6 


2.90 x 10~ 5 


4.57 x 10~ 2 


9.88 x 1(T 4 


1.19 x 10" 1 


10- 7 


3.46 x 10~ 6 


6.90 x 10~ 3 


6.87 x 1(T 5 


2.00 x 10~ 2 


1(T 8 


3.85 x 1(T 7 


7.92 x 10" 4 


6.71 x 1(T 6 


2.26 x 10~ 3 


10 -9 


3.63 x 1(T 8 


1.01 x 10" 4 


9.16 x 1(T 7 


2.60 x 10~ 4 


io- 10 


2.71 x 1(T 9 


1.28 x 10~ 5 


6.35 x 1(T 8 


3.10 x 10~ 5 



Table 3: r — 6; sin8(V r ,V r ) and sin6(W r , V r ) as s varies 



4. Backward error 

Instead of seeking bounds on how much an inexactly computed reduced 
model differs from an exactly computed counterpart, one may view an inex- 
actly computed reduced order model as an exactly computed reduced order 
model of a perturbed full order system. That is, we wish to find a full order 
system 

•K(s) = e(s)X(s)- 1 S(s) (36) 

so that the inexactly computed reduced model for 9£(s) = Q(s)DC(s)~ 1 'B(s) 
would be an exactly computed interpolatory reduced model for Ji(s). Given 
left and right tangential interpolation data as in (J9J) and (10) that has con- 
tributed toward producing the inexactly computed interpolatory reduced 
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model U£ r (s), find IK(s) as in (36) so that 

cfiK(jAi) = cJO-Cr(fii) for i — 1, . . 
fK(<Tj)bj = IK r (cr.,)b.,- for j = 1, 



r, and 
. . , r. 



and so that IK r could have been computed from the perturbed system IK from 
the given tangential interpolation data via an exact computation. Specifi- 
cally, given computed (inexact) projecting bases 



V r = [vi, 



w 



w 



as in (15), and a resulting (inexact) reduced order coprime realization 

•k r (s) = e r (s)3c r (s)" 1 s r ( s ), 

find a full-order system IK(s) = G(s)JC(s)~ 1 'B(s) so that left and right inter- 
polation conditions hold: 

cfe(ja) = wfXQii) for i 
^i)vj = &(<rj)bj for j 

and so that 

W^K(s)V r , B r (s) = W r T S(s), 



(37) 
(38) 



UC(s) 



and ejs) = C(s)V r 



(39) 



There (typically) will be an infinite number of possible systems, IK, that 
are consistent with the computed reduced system IK r in this sense - - we 
are interested in those that are close to the original system IK with respect 
to a convenient system norm such as "H^ or % 2 - in order to proceed, it is 
convenient to restrict the class of backwardly compatible systems, IK. We 
consider those that have realizations that are constant perturbations from 
the corresponding original system factors: 

X(s) = X(s) + F, S(s) = S(s) + E, and C(s) = C(s) + G. (40) 



where E, F, and G are constant matrices. The conditions (37), (38), and 



(39) impose constraints on E, F, and G. Indeed, (37) and (38) imply that 

T 



wf F + £ 
Fv, + ry, 



cfE for i = 1, 
Gbj for j — 1, 



r, and 
, r. 
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(39) implies that 



W T T FV r = 0, WjG = 0, and EV r 



0. 



Taken together, we find that backward perturbations of the form (40) can 
exist only if 



£jy r = for i 



and Wjrj,: = for j — 1, 



(41) 



Thus, we find constraints on the inexact interpolation residuals $, t and rj in 
order for a backwardly compatible system of the form (40) to exist. More 



complicated perturbation classes than (40) may be considered that would 



allow us to remove the conditions (41), of course, but instead we choose to 
focus on a computational framework that guarantees (41). The Biconju- 



gate Gradient Algorithm (BiCG) will be an example of an iterative solution 
strategy that fits this framework [TJ 0]; others can be constructed without 
difficulty, although many standard strategies such as GMRES, do not fit this 
framework. 



4-1. The Petrov-Galerkin Framework for Inexact Solves 



We have observed above that (41) is necessary for there to be a well- 



defined backward error of the form (40) to exist. The simplest framework 



within which one may generate reduced order models that are guaranteed 
to satisfy this condition involves a Petrov-Galerkin formalism for produc- 



ing approximate solutions to (11) and (12). For simplicity, we restrict our 
discussion to the case that ^ = ai (identical left and right interpolation 
points). 

Let Vn and Qn be iV-dimensional subspaces of C n satisfying a nondegen- 
eracy condition: (JC^ct^Vn) H Qn = {0} for all shifts, <Tj to be considered. 
The Petrov-Galerkin framework for generating approximate solutions to the 



interpolation conditions (11) and (12) proceeds as follows 



Find v,- e Vn so that DC(cr,)v, 



Sfo)bi 



Qn and 



find Wj e Q N so that 3C(<jj) T Wj — C(<7j) T Cj _L Vn 



(42) 



Computed quantities generated within a Petrov-Galerkin framework will be 
denoted with a "tilde" to distinguish them from earlier "hat" quantities where 
no structure was assumed in the inexact solves. The following theorem asserts 
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that if a reduced order model is computed within a Petrov-Galerkin frame- 



work (42), then one can obtain a structured backward error that throws the 
effect of inexact solves back onto a perturbation on the original dynamical 

system. 

Theorem 4.1. Given a full order model IK(s) = C(s)3C(s) _1 2J(s), inter- 
polation points {ctjYj = i, and tangent directions {bj}[ =1 and {cj}[ =1 , let the 
inexact solutions Vj for X(crj)~ le B(aj)hj and w^for JC(aj)^ T G(aj) T Cj be 



obtained in a Petrov-Galerkin framework as in (42). Let\ r and W r denote 
the corresponding inexact interpolatory bases; i.e. 

V r = [v 1 , ••-, v r ] and W r = [ w 1; ■■-, w r ] . (43) 

Define residuals 

Vj = OCicr^vj - SfoOb,- and £,. = Xfafwj - efafcj, 

residual matrices 

Rb= [Vi, V2, •••,»7r]» R c = [£i, £ 2 , •••>&.]> ( 44 ) 

and the rank 2r matrix 

F 2r = RbtWfCg^Wj + VrCW^Vr)- 1 !^. (45) 

Let tK r (s) = G r (s)0C r (s)~ lt B r (s) denote the computed inexact reduced model 
via the Petrov-Galerkin process where 

K r (s) = WjX(s)V r , S r (s) = W^B(s), and G r (s) = C(s)V r . (46) 

Then, J€ r (s) exactly tangentially interpolates the perturbed full- order model 

5i{s) = e(s){X(s) + F 2r )- 1 3(s), (47) 

at each c^: 

and cf !K (oi)bj = cf IK r (crj)bj for each i — 1, . . . , r. 
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Proof: The computed model, IK r (s), will (exactly) tangentially interpolate 
a perturbed model Oi(s) = G(s)(X(s) + F) -1 3(s) provided the following 
interpolation conditions hold: 

(X(<Ti) + F) v, = 3(ai)bi and w[ (X(ffi) + F) = c^Cfo) for i = 1, . . . , r. 

Equivalently, these can be interpreted as conditions on the perturbation F. 
Rewriting this using notation defined above, F must satisfy 

FV r = R b and W^F = R£. (48) 

The Petrov-Galerkin framework guarantees W^Rb = and R^V r = 0. 



Substitution of F 2r from (45) into (48) verifies that F 2r is a perturbation to 



X(s) for which the computed (inexact) vectors become (exact) interpolation 
vectors. 

Note that since W r T F 2r V r = 0, 

X r (s) = WjX{s)V r = Wj(X(s) + F 2r )V r . 

Consequently, the reduced model Ji r (s) obtained by inexact solves in (461) 
is what one would have obtained by exact interpolatory model reduction of 
W(s). □ 

Theorem 4.2. Assume the hypotheses of Theorem Ujl and that W^V r is 
nonsingular. Define an oblique projector, $ r = V r (WjV r ) _1 W^. The 
backward perturbation F 2r given in Theorem [O satisfies 



||F 2r || F < Vr||* r || • (maxP4^min(V r D)- 1 + max|l4^min(W r D) l ) 
V l ll v *ll i \\ w i\\ J 

where <^ m ; n denotes the smallest singular value and \\M-\\f = ytrace(M^M) 
denotes the Frobenius norm of a matrix, M. 

Proof: Note that 

||F 2r || F < ||R b (W r T V r )- 1 W r T || F + HV^WjV,.)- 1 ^^. 
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Let V r have an orthogonal factorization as V r = Q^L^ with Q*Qt, = I. 
Then 



|R b (W;!V r )- 1 W;{|| F 



= 


IRbL^^^W^V,.)- 1 ^!^ 


< 


|R b L- 1 || i ,.||L t; (WjV r )- 1 Wj 


< 


RbL~ \\f ■ \\&r\\ 


< 


IRbD^L.D,)- 1 !^-!!^!! 


< 


IRbD^i^-iKL^)- 1 !!-!!^!! 



where we have introduced a diagonal scaling matrix 

D v = diag(l/||vi||, l/||v 2 ||, ..., l/||v r ||). 

~ i~ ll^'ll 

Easily one sees RbDJIp < y/rmax \ . For the remaining term, note that 

» v,- 



(L„D„) 



V r D,,x| 



mm 



Mnin I * r*-* 



A similar bound for || V r (W;TVV) _1 Rc ||f is produced by an analogous process, 
which leads then to the final estimate for ||F 2r ||^. □ 

Note that the perturbation F 2r is completely determined by accessible, 
computed quantities. Hence, one can use F 2r to determine how accurately 
one must solve the underlying linear systems in order to assure system fidelity 
of a given order. 

Theorem 4.3. // ||F 2r || < l/||3C(s)- 1 || Woo then 



\<K{s)-K{s)\\ H2 < 



lOCfs)- 1 ! 



He 



"2r 



'2r 



PROOF: The system-wise backward error associated with inexact solves 
may be written as 



F 2r y 1 'B(s) 



M(s) - M(s) = e(s)X(s)- 1 'B(s) - e(s) (X(s) 
= e(s)K(s)- 1 F 2r (X(s)+F 2r )- 

= e(s)3C(s)- 1 F 2r (l + 3C(sy 1 F 2r y i JC(s)- 1 'B(s) 



F2r) _1 3(s) 

-l 
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Define M(s) = F 2r 

m(s)-ns)\\ 2 H2 = 



I + DC(s) 1 F 2r ) and observe that 



2tt 



< 



2tt 

1 

2^ 



leM^M" 1 !!! 



oo 

oo 



\e{iu)'K{iuY 1 ^i{iuj)'K{iujY 1 'B{iu)\\ 2 F du 
M(zw) f ■ WX^iujY^^W 2 du 

< I — / liei/^lJKl/^r'lli-r/u;! • max||MM f • max ||3C(«^) _1 S(«07) || 2 

< \\G(s)0C(s)-Tn 2 ■ \\X(*)- l *(*)\\k. ■ II^^HL- 

To estimate H^Vl^s)!!?^,, a rearrangement of the definition of M(s) provides 

M{s) = (l-M(s)0C{s)- 1 )F 2r . 
So we have immediately, 

\\M(s)\\ Hoo =max\\M(tuj)\\<max\\I-M(iuj)OC(iuj)- 1 \\ ■ ||F 2r || 

< tl + max\\M(iuj)3C(iujy l \\ J ||F 2r || 
<(l + \\M(s)\\ Hoo HXCa)- 1 !!^) ||F 2r || 
Since IIDC(s)^ 1 H^ ||F 2r .|| < 1, this last expression can be rearranged to obtain 

IITP II 



||M( S ) 



•Ho 



< 



|3Cfs) 



Ho 



|F 2r j 



which implies the conclusion. □ 

By combining Theorem |3.3 with Theorem 3^ or combining Theorem 4J2 



with Theorem 4J3 we approach our goal of connecting quantities that we 
have control over, such as the termination threshold, e, to relevant system 
theoretic errors, \\Ji r — 3€ r || and ||J£ — 3f ||, which are quantities we would 
like to control. 

One may use these expressions as a basis to devise and investigate differ- 
ent, effective stopping criteria in large-scale numerical settings. For example, 



while e appears explicitly in Theorem 3.2 in a way that suggests its use as a 



relative residual norm threshold; while Theorem 4.2 suggests a scaling of the 
residual norm by the norm of the solution vector as another possible stopping 
criterion. These and related ideas are the focus of on-going work. 
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4-2. Quantities of interest in derived bounds 



By combining Theorem 4.2 with Theorem 4.3, one observes that pertur- 
bation effects of the inexact solves on the system theoretical (model reduction 
related) measures critically depend on the four quantities: The norm of the 
oblique projector <fr r = V r (WjV r ) _1 W^ of the underlying model reduction 
problem, reciprocals of the minimum singular values of the scaled primitive 
bases V r D and W r D; and the stopping criterion e for the inexact solves, 

11^7'H ll^'ll 

(which affects max l and max ' l 



JVi|| i ||Wi|| 

The e term is associated directly with inexact solves and is under the 
control of the user. The remaining quantities SininCVrD) -1 , ^ min (W r D) _1 
and ||3> r ||, depend largely on the selection of interpolation points {cTi} and 
tangent directions, but the influence of interpolation data on the magnitude 
of these quantities is difficult to anticipate. 

In this section, we will investigate experimentally the effects of the in- 
terpolation point selection on the three quantities of interest, Wn(V r D) -1 , 
<Jmin(W r D) -1 and || $7-11, appearing in the derived bounds. These quantities 
are continuous with respect to the primitive basis vectors, {vi, • • • , vy} and 
{wi, • • • , w r } in neighborhoods where W^TV r is nonsingular (i.e., where the 
projector <& r is well defined). Thus it will be sufficient to examine how the 
magnitudes of the quantities of interest depend on interpolation data presum- 
ing that the necessary linear solves are done exactly; for modest convergence 
thresholds, the effect of inexact solves on these magnitudes is secondary to 
the effect of interpolation point location. 

For our numerical study, we use the International Space Station 12A 
Module as the full-order model. The model has order n = 1412. We exam- 
ine a single-input single-output subsystem, ^K(s), reducing the order from 
1412 to order r with r varying from 2 to 70 in increments of two. For 
each reduced order, we chose 2000 random shift selections and computed 
<smin(V r D) -1 , <j min (W r D) _1 and ||<fr r ||- For each r, r/2 shifts were sampled 
from a uniform distribution on a rectangular region in the positive half- 

min A |Re(A)|<Re(2;)<max A |Re(A)| \ 

i, , i ' II /am >, where the max and 

|lm(z)| < max A |lm(A)| J 

min are taken over all the poles of the system. The remaining r/2 shifts were 

taken to be the complex conjugates of this random sample, so as to produce a 

shift configuration that was closed under conjugation. Additionally for each 

r, we applied model reduction using the "H 2 -optimal interpolation points gen- 
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erated by the method of [TB|. Then, for each r, out of the 2000 randomly 
generated shift selections, we counted the number of cases where the ran- 
dom shift selection yielded smaller values of ^ m i n (V r D) _1 , g m i n (W r D) _1 and 
||<fr, r ||. The results are shown in Figure [2] Figure [2]-(a) and -(b) show that for 
most of the cases, the %2-optimal interpolation points yield smaller values 
for Win(V r D) _1 , <r m i n (W r .D)~ 1 . Indeed, for r > 48, the %2-optimal points 
produced smaller values in more than 99% of the cases. Also, for the last 
three cases: r = 66, r = 68, and r = 70, the "^-optimal interpolation points 
always yielded smaller quantities. The results are even more dramatic for the 
projector norm, which is important in scaling the perturbation effects caused 



by inexact solves, see Theorem 4.2 Out of 70, 000 cases (2000 selections for 
each r value), the %2-optimal interpolation point selection produced smaller 
condition numbers in all except 7 instances: 5 instances for r = 2, and 2 
instances for r = 8. These numerical results illustrate that ^-optimal in- 
terpolation points can be expected to yield smaller values for ^ min (V r .D) _1 , 
g min (W r .D) _1 and ||«fr r ||, and hence should produce reduced order models 
that are more robust with respect to perturbations. 

Figure [2] also shows that for r = 14, 48% of the randomly selected shifts 
yielded smaller values of ^mmCVrD)" 1 . However, when we inspected the 2000 
randomly selected shift sets for r = 14 in more detail, we observed some inter- 
esting additional features. We computed the three quantities ^ min (V r D) _1 , 
Wm(\V r D) -1 a nd ||«fr r || for each of the 2000 randomly selected shift sets, and 
compared them with the corresponding value derived from an ^-optimal 
shift selection. The results are shown in Figure [3j The top plot shows 
Wn(V r D)/^ min (V° pt D) where V° pt stands for the primitive interpolatory 
basis for the T^-optimal points. The bigger this ratio, the better the random 
shift selection. Even though for 48% of the cases, the random selection was 
better, the highest this ratio becomes is 2.20, i.e., the random shifts were 
never much better than a factor of 2 better than what "^-optimal shifts 
provided. For the remaining 52% of the cases, the randomly selected shifts 
were worse, and often worse by a factor of 100 or more. The situation for 
W r is shown in the middle plot. Once more, the situation is much more 
drastically in the favor of the ^-optimal interpolation points when the pro- 
jector norm is inspected; the bottom plot in Figure [3] which depicts the ratio 
||$ r ||/||$° pt || where $° pt denotes the projector for the T^-optimal points. 
As illustrated in Figure [2j there are no random shift cases yielding a smaller 
projector norm. Furthermore, in many cases the projector norm for the ran- 
dom shift selection is almost 4 order of magnitudes higher than that of the 
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Figure 2: Comparison of ? m j n (V r D) 1 , <j min (W r D) x and ||# r || for random shift selections 
relative to values for T^-optimal shifts 



T^-optimal points. Indeed, on average the projector norm for the random 
points is 8.19 x 10 1 times higher. These numbers change more in the favor of 
the T^-optimal points as r increases. For example, for r = 50, while the ratio 
Wn(V r D)/^ min (V° pt D) becomes only as high as 1.48, it becomes as low as 
2.89 x 10~ 4 for some random selections; Also, the ratio 1 can reach as high as 
2.91 x 10 5 . For r = 70, \\*& r \\ for random selection is 1.73 x 10 2 times higher 
than ||*fr° pt || on average. The three quantities we have been investigating 
appear to be extremely well conditioned for - H 2 -optimal interpolation points. 



Even for r = 70, both Win(V° pt D) 
and ||$° pt || is smaller than 7. 



Wn(W,° pt D) x remain smaller than 10 
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Figure 3: Detailed comparison for r = 14 



5. Inexact Solves in Optimal Interpolatory Approximation 

The quality of the reduced-order model in interpolatory model reduction 
clearly depends on the selection of interpolation points and tangent direc- 
tions. Until recently, this selection process was mostly ad hoc, and this factor 
had been the principal disadvantage of interpolatory model reduction. For 
systems in standard first-order state-space form, Gugercin et al. [16] have 
produced that an T^-optimal interpolation point / tangent direction selec- 
tion strategy and proposed an Iterative Rational Krylov Algorithm (IRKA) 
to generate interpolatory reduced-order models that are (locally) optimal 
with respect to the H2 norm. (An "^-optimal interpolation point selection 
strategy is still unknown for the general coprime factorization framework.) 
In this section, we investigate the behavior of inexact solves within the %2- 
optimal interpolatory approximation setting, specifically examining the be- 
havior when inexact solves are employed in IRKA. In the rest of this section, 
we briefly review the optimal H2 approximation problem and the method of 
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We then show how inexact solves can be employed effectively in this 
setting and discuss observed effects on optimality of the final reduced model. 
Our discussion focuses on systems in first-order descriptor form: 

!K(s) = C(sE - A) _1 B (49) 

where E, A e R nxn , B G R nxm and C G R pxn . 

5.1. Optimal H2 approximation problem 



Given the full-order system as in (49), the goal of the optimal [K 2 model 



reduction problem is to find a reduced-order model J£ r (s) that minimizes the 

"H 2 error; i.e. 

IIJC-JCJL = min IIIK-GUL . (50) 

11 rUH ' 2 G r stable " rUH2 v ; 

dim(G r )=r 

Many researchers have worked on this problem. These efforts can be 
grouped into two categories: Lyapunov-based optimal "H 2 methods such as 
[3"Tl [251 [T71 [THl EHl E2] ; and interpolation-based optimal "H 2 methods such as 
[23l[l6l[T5l[29l[ini[Tll2ni[5l[7]. Here, we will focus on the interpolation- 
based approach. However we note that Gugercin et al. [16] has shown that 
these two frameworks are theoretically equivalent; hence motivating the use 
of interpolatory approaches to optimal % 2 approximation since they are nu- 
merically superior to the Lyapunov-based approaches. 



Since the optimization problem (50) is nonconvex, obtaining a global 



minimizer is a hard task and can be intractable. The usual approach is 
to find reduced order models that satisfy first-order necessary optimality 
conditions. Meier and Luenberger [23] introduced interpolation-based % 2 - 
optimality conditions for SISO systems. Analogous "H 2 -optimality conditions 
for MIMO systems have recently been developed by [HJl [ITJJ [2H] which in turn 
have led to analogous algorithms for the MIMO case; see [T6J, [TO] for more 
details. 

Theorem 5.1. Given "K(s) = C(sE - A) _1 B, let % r (s) = E[=i ^T c i h I 
be the best r th order approximation of IK with respect to the % 2 norm. Then 

(a) K(-\ k )b k = J€ r (-A fc )b fc , (b) c£K(-A fc ) = cJlK r (-A fe ), (51) 
and (c) cl'K' (-\ k )h k = cl^K' r (-\ k )h k for k = 1, 2, 



j, ..., 



r. 
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5.1.1. An algorithm for interpolatory optimal rl<i model reduction 



Theorem 5.1 reveals that any rli optimal reduced-order model ZK r (s) is 
a bi-tangential Hermite interpolant to !K(s) at mirror images of the reduced- 
order poles. However, since the interpolation points and the tangent direc- 
tions (and consequently, V, r and W r ), depend on the final reduced-model to 
be computed, they are not known a priori. The Iterative Rational Krylov 
Algorithm (IRKA) of [16] resolves this problem by iteratively correcting the 
interpolation points and the directions as outlined in Algorithm [TJ The 
reduced-order order poles are reflected across the imaginary axis to become 
the next set of interpolation points; the tangent directions are corrected using 
residue directions from the current reduced model. Upon convergence, the 
resulting interpolatory reduced-order model satisfies the necessary conditions 



of Theorem 5.1 For further details on IRKA, see 116 



Algorithm 1. IRKA for MIMO V.2 Optimal Tangential Interpolation 

1. Make an initial r-fold shift selection: {a±, . . . , oy} and initial tangent direc- 
tions h>i , . . . , b r and 61 , . . . , c r . 



2. V 



r 



r 



(<tiE - A)- 1 Bbi ••• (<r P E - A) _1 Bbr 
W r = [((Ti E - A^-^C! • • • (oy E - A T )- 1 C T ci 

3. while (not converged) 

(a) A r = WjAV r; E r = W^EV r; B r = W^B, and C r = CV r 

(b) Compute Y* A r X = diag(Xi) and Y*E. r X = I r where Y* and X are 
the left and right eigenvector matrices for AE r — A r . 

(c) o~i < — — Aj(A r ,E r ) for i = l,...,r, b* < — e^Y*B r and Cj < — 

(d) V r = (aiE-A^Bbi ••• (a r E - A)- x Bb r 

(e) W r = [(en E - A T y 1 C T c 1 • ■ ■ (oy E - A T )~ 1 C T c 1 ] . 

4. A r = WjAV r , E r = W^EV,., B r = W^B, C r = CV r 

5.2. Inexact Iterative Rational Krylov Algorithm (InxIRKA) 

For large system order, one may see from Algorithm [TJ that the main cost 
of IRKA will generally be solving 2r large linear systems at each step. If the 
IRKA iteration converges in k steps, a total of 2rk linear systems will need 
to be solved. In settings where system dimension reaches into the millions, 
iterative linear system solvers become necessary and inexact linear system 
solves must be incorporated into IRKA. We refer to the modified algorithm 
as the Inexact Iterative Rational Krylov Algorithm (InxIRKA) and describe it 
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in Algorithm [2] below. We employ the Petrov-Galerkin framework for the 
inexact solves. In Algorithm [2j the function Fpc in 

[v,, w,] = F PG (A, E, B, a u b h c h v (0) , w( 0) , e) 

denotes an inexact solve using a Petrov-Galerkin framework to approximately 
solve the linear systems (<TjE — A)vj = Bbj and (ctjE — A) T Wj = C T Cj with 
initial guesses v^ ) and w^ \ respectively, and a relative residual termination 
tolerance e, i.e., at the end, 

II (<jiE - A) v 4 - Bbi II ,|| (ct,E - A) T w, - C T Cl || 

"<e and ^ „' , -<£ 



|Bb;|| - ||C T C; 



Algorithm 2. InxIRKA for MIMO H2 Optimal Tangential Interpolation 

1. Make an initial r-fold shift selection: {a"i, . . . , oy} and initial tangent direc- 
tions bi , . . . , b r and c±, . . . ,c r . 

2. [vj, Wj] = F PG (A, E, B, en, b i; Cj, 0^0, g) for i = 1, . . . , r 

3. V r = [vi, v 2 , ...,v r ] and W r = [ wi, w 2 , ...,w r ]. 

4. iu/ji/e fnoi converged) 

(a) A r = W^AV r , E r = W^EV r , B r = W^B, and C r = CV r 

(b) Compute Y*A r X = diag(Xi) and Y*E r X = I r where Y* and X are 
the left and right eigenvector matrices of AE r — A r . 

(c) Oi < Aj(A r , E r ) for i = l,...,r, b* < — ej Y*B r and 

(d) [vj, Wj] = Fp G (A, E, B, ctj, bj, c^Vj, Wj, e) for i = 1, . . . ,r 

(e) V r = [ vi, v 2 , . . . , vy ] and W r = [ wi, w 2 , . . . , w r ]. 

5. A r = W^AV r , E r = Wj£V r , B r = W^B, and C r = CV r 

As discussed and illustrated in [HI [3], in most cases IRKA converges 
rapidly; that is, the interpolation points and directions at the k th step of 
IRKA stagnate rapidly with respect to k. Let a\ and bj denote the i th in- 
terpolation point and right-tangential direction, respectively, at the k th step. 
Then we expect that as k increases, the solution vj of the linear system 
(crj E — A)v( fc ) = Bbj from the k th step approaches to the solution vj 
of the linear system (cxf +1) E- A)v( fc+1) = Bbf +1) at the (k + l) st step. This 
is precisely the reason that in Step 4.(<f) of Algorithm K2J we use vj as an 
initial guess in solving (crf +1) E - A) V ( fe+1 ) = Bbf +1) at the (k + l) st . We 
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expect that this initialization strategy will speed-up the convergence of the 
iterative solves. 

The development of effective stopping criteria based rationally on sys- 
tem theoretic error measures as we have introduced them here is the focus 
of on-going work. Similar approaches toward the design of effective precon- 
ditioning techniques and reuse of preconditioners tailored for interpolatory 
model reduction and especially for optimal H2 approximation are also under 
investigation. 

5. 3. Effect of Inexact Solves in the InxIRKA Setting 

The first question to answer in InxIRKA is whether a statement can be 
made about the optimality as in the exact IRKA case. Employing the Petrov- 
Galerkin framework makes this possible: 

Corollary 5.1. Let !K r (s) be obtained by Algorithm 2, Then 3£ r (s) satisfies 
the necessary conditions for optimal %2 approximation of a near-by full- order 
model 9€(s) = C(sE — (A + F2 r ))~ 1 B where F2,. is the rank-2r perturbation 



matrix defined in (45) 



Corollary |5.1| shows that with the help of the underlying Petrov-Galerkin 
framework, we state that the final reduced model of InxIRKA is an optimal 
%2 approximation to a nearby full-order model. 



As we discussed in Section |42| for a good selection of interpolation points, 
interpolatory model reduction is expected to be robust with respect to per- 
turbations due to inexact solves. Hence, if one feeds the optimal interpolation 
points from IRKA into an inexact interpolation framework, we expect that the 
resulting reduced model will be close to the optimal reduced model of IRKA. 
However, the optimal interpolation points are not known initially and In- 
xIRKA will be initiated with a nonoptimal initial shift selection. If the initial 
interpolation points and directions are poorly selected, at the early stages 
of the iteration, perturbations due to inexact solves might be magnified by 
this poor selection. One can avoid this scenario by using a small termination 
threshold e in the early steps of InxIRKA, and then gradually increase e as 
the iteration starts to converge. However, we note that in our numerical ex- 
periments using random initialization strategies, InxIRKA performed robustly 
and yielded high fidelity reduced models that are also close to the true opti- 



mal reduced model. This is illustrated in [5A_ below. Effective initialization 
strategies are discussed in [TB] as well. 
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5-4- Numerical results for InxIRKA 

Here we illustrate the usage of inexact solves in the optimal H2 approx- 
imation setting by comparing IRKA with InxIRKA. We use the example of 



£ 3.3 but with a finer discretization leading to a state-space dimension of 
n = 20209. We focus on a MIMO version using 2-inputs and 2-outputs. 

We reduce the order to r = 6 using both IRKA and InxIRKA. In InxIRKA, the 
dual linear systems are solved in a Petrov-Galerkin framework using BiCG 
[I] where we use three different values for the relative residual termination 
threshold of e: 10~ 5 , 10~ 3 , and 10 _1 . In all cases, the behavior of InxIRKA is 
virtually indistinguishable from that of IRKA. Starting with the same initial 
conditions, both IRKA and InxIRKA converge within 10 iteration steps in all 5 
cases. The evolution of the T-L 2 errors ||IK — 3~C r ||% 2 and ||IK— JK r ||% during 
the course of IRKA and InxIRKA, respectively, are depicted in the top plot of 
Figure |4j The figure shows that InxIRKA behavior is almost an exact replica 
of that of IRKA. The deviation from the exact IRKA is noticeable in the graph 
only for e = 10 _1 . To illustrate how much ^K r deviates from !K r as IRKA and 
InxIRKA evolve, we show the progress of \\0€ r — 3~C r ||« 2 i n the bottom plot 
of Figure |4j For this example, we initialized both IRKA and InxIRKA with 
an initial reduced-order model (as opposed to specifying initial interpolation 
points and tangent directions). Thus, IK r = c H. r initially and no linear solvers 
are involved in the first (k = 0) step. One could expect that perturbation 
errors due to inexact solves might accumulate over the course of the InxIRKA 
iteration, but this does not appear to be the case as this figure illustrates. 
The magnitude of ||!K r — IK r ||^ 2 remains relatively constant throughout the 
iteration at a magnitude proportional to the termination criterion. 

The resulting "H 2 and "Hoc, model reduction errors, ||U-C — 3~Cr ||w 2 and 
||J€ — 3^r||woo (with !K r obtained from IRKA), versus ||IK — JK r ||« 2 and 
\\0i — ^K-tWhoo (with !K r obtained from InxIRKA) are given as e varies in 
Table [4] below. The row corresponding to e = represents the errors due 
to exact IRKA. These numbers demonstrate that employing inexact solves in 
InxIRKA does not degrade the model reduction performance. We also mea- 
sure the difference between !K r and "Ji r in both H2 and Hoc, norms as e 
varies. These results are tabulated in Table [5] Note that while ||IK — 3~C r ||^ 2 
and ||3~C — ^Crllw^ are respectively C(10 -4 ) and O(10~ 2 ), the contributions 
attributable to !K r — !K r are much smaller in magnitude and do not al- 
ter the resulting (optimal) model reduction performance in any significant 
way. If one were to convert the perturbation errors in Table [5] to relative 
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Figure 4: Evolution of the "H 2 error during IRKA and InxIRKA 

error (as opposed to the displayed absolute error), both \\5C r — J-C r \\ n2 and 
\\^K r — ^Krllwoo starts at (9(10~ 6 ) for e = 1CT 5 , and increases linearly by one 
order as e increases by the same amount. 

We finally list, in Table [6j the final exact and inexact optimal interpo- 
lation points due to IRKA, and InxIRKA for e = 10~ 3 and e = 1CT 1 : Not 
surprizingly, the resulting interpolation points are very close to each other 
(though not the same). This can be viewed as another illustration of the fact 
that 3£ r is an %2 optimal approximation to a nearby full-order system. 

As discussed above, in the implementation of InxIRKA, we used the solu- 
tion vectors from the previous step as the initial guess for the linear system in 
the next step taking advantage of the convergence in the interpolation points 
and tangent directions. To illustrate the effectiveness of this simple approach, 
throughout InxIRKA we monitor the number of BiCG steps required to solve 
each linear system. We illustrate the behavior only for one of the interpo- 
lation points. We choose the interpolation points closest to the imaginary 
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3.708415753 x 10" 4 


1.084442854 x 10" 2 


io- 5 


3.708415754 x 10~ 4 


1.084425703 x 10~ 2 


10" 4 


3.708415778 x 10" 4 


1.084282001 x 10~ 2 


10- 3 


3.708418102 x 10~ 4 


1.082437228 x 10~ 2 


10- 2 


3.708621743 x 10~ 4 


1.064836300 x 10~ 2 
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3.716780975 x 10" 4 


1.055441476 x 10~ 2 



Table 4: Evolution of the model reduction errors as e varies 
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IIT-f r \-C II 
|| j\, r */v r ||% 2 


IIT-T C VC II 
\\JX r — J^rWHoo 


10~ 5 


5.1921 x 10~ 9 


2.7776 x IO" 7 


10- 4 


5.7156 x 10~ 8 


2.4611 x 10~ 6 


10~ 3 


6.3982 x 10~ 7 


2.1043 x 10~ 5 


10- 2 


5.9277 x 10~ 6 


2.0910 x 10~ 4 


10- 1 


2.2056 x 10~ 5 


2.9228 x 10~ 3 



Table 5: Evolution of the perturbation error as e varies 

axis since these produce the hardest linear systems to solve and invariably 
contribute most to the cost of inexact solves. Figure [5] depicts the the num- 
ber of BiCG steps required as InxIRKA proceeds for these interpolation points 
using three different stopping criteria e = 10~ 5 , e = 10~ 3 and e = 10 _1 . The 
figure clearly illustrates that re-using the solutions from the previous steps 
works very effectively in reducing the overall cost of the BiCG. The number 
of BiCG steps goes from 1200 down to 200 in 3 to 4 steps. 



Evolution of BiCG step effort 
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Figure 5: Evolution of BiCG effort during InxIRKA for shift closest to the imaginary axis 
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Table 6: Optimal interpolations points as e varies 



6. Structure-preserving interpolation for descriptor systems 

The backward error analysis of ^4] has been presented for the transfer 
functions in the generalized coprime factorization form as in pi). In this 
section, we show that stronger conclusions on the structure of the reduced 
system can be drawn in the case the system has a realization as a descriptor 
system, that is, 

M{s) = C(sE - A^B (52) 

where E, A G R nxn , B G E nxm , and C G R pxn are constant matrices. In this 
case, for the interpolation points {oj} r . =v and the tangent directions {bj}^ =1 
and {cj} r . v the associated primitive interpolatory bases V. r and W r can 
be obtained from (13) and (14) using 3C(s) = sE — A, S(s) = B (constant 
matrix) and G(s 



C (constant matrix). Then, the resulting reduced-order 



model is given by 






V_; r yShj r 



X\.rf I Jjn- 



where 

E r 



WrEV r , A r = WjAV r , B r = WrB, and C r = CV r 



(53) 



(54) 



Let the set S = {<Ji, bj, Cj} denote given tangential interpolation data. De- 



fine the matrices LpK,«S] G C rxr and M[JK,5] G C 
dynamical system !K(s) and interpolation data S: 



corresponding to the 



CMK,s]) y 



cJ{-K{a l )-'K{a ] ))h J 



<?% - o~j 



cTH'laJbi 



if i^j 



if 



(55) 
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cf (di'Ktoj) - <TjK{<jj)) bj 



(Mpi,S]) 



Oi 



(Ji 



'■j 



of [sn^'is^M 



if i^j 



if 



(56) 



L[!K, <S] is the Loewner matrix associated with the interpolation data S and 
the dynamical system IK(s), MpK, S] is the shifted Loewner matrix asso- 
ciated with the interpolation data S and the system sfK(s), see [31 122] ■ 
The next theorem presents a canonical structure for the exact interpolatory 



reduced-order model (53)-(54). 

Theorem 6.1. l2~2^ Given a full-order model !K(s) = C(sE — A) _1 B and 
tangential interpolation data S = {o~i, bj,Cj}, i/ien t/ie reduced-order quanti- 



ties in (54) satisfy 



E r = — 1L[!J~C, oj, 
A r = -MpK,«S], 

C r = [ JK((Ti)bi, 



B, 



c^(a r ) 
IK(<7 r )b r ]. 



(57) 



£.i. T/ie Petrov-Galerkin framework and structure preservation 



Theorem 6.1 presents a canonical form for the exact bitangential Hermite 



interpolant in the case of standard state-space model. Next we show that if a 
Petrov-Galerkin framework is employed in the solution of the linear systems, 
the inexact reduced-model will have exactly the same form as the exact one. 
The result is a direct consequence of Theorems 4^ and |6.1| 



Corollary 6.1. Given the standard full- order model !K(s) = C(sE — A) X B 
together with the the interpolation data S = {<7j,bj,Cj}, let the inexact solu- 



tions Vj for (ctjE — A) Bbj and w ? - for (o~jE — A 



~ T C T Cl 



be obtained in 



a Petrov-Galerkin framework as in (42). Let V r and W r denote the corre- 



sponding inexact Krylov bases as in (43). Define the residuals 



Bb, 



and 



rfj = (ctjH; - AjVj - ttbj and £^ = (cr,E - A) w,- - C c,. 
Let the residual matrices Rb and R c , and the rank 2r matrix F2,. be as defined 



in (44) an d (45), respectively. Then, the inexact interpolatory reduced-order 

5< r (s) = C r (sE r - A r ) _1 B r (58) 



model 
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is an exact Hermite bitangential interpolant for the perturbed full- order model 



M(s) = C(sE - (A + Far))"^. 
Morever, the reduced-order quantities satisfy 



E r 



-LpK,S], 
■MpK,S], 



B, 



_ c^'K{a r ) _ 



C 



M((n)bi, 



(59) 



^(a r )b r ].(60) 



where LpK, <S] and M[IK, «S] are i/ie Loewner matrices associated with the 
dynamical systems Ji(s) and sfK(s) respectively, and the interpolation data 



S as defined in (55) and (56) 



Corollary 6T reveals that the inexact reduced-order model quantities have 
exactly the same structure as their exact counterparts. The interpolation 
data S is the same in both cases; the only difference is that 3€(s) is replaced 
by U£(s) in the construction that yields the Loewner-matrix structure. The 
preservation of this structure is independent of the accuracy to which the 
linear systems are solved. In the case where E = I, the structure of the exact 
and inexact reduced-models becomes even simpler: 

Corollary 6.2. Assume the hypotheses of Theorem 6A_ with E = I. 
the exact interpolant tK r (s) = C r (sl r — A r ) _1 B, r satisfies 

A r = £-QB, B r = Q, and C r = [ H(<Ti)bi, . . . , H(a r )b r 

where 



Then 



(61) 



Q = (W^V r ) _1 W^B, £ = diag(cr 1 , ...,a r ) and B = [b 



i- 



,b r ]. (62) 



Assume the hypothses of Corollary 6.1 with E 
polant Jir(s) = C r (sl r — A r ) _1 B r . satisfies 



A r = £ 
where 



QB, B r = Q, and C r 



I. Then, the inexact inter- 
[9C(<ri)bi, ..., J€{a r )b r ] (63) 



Q = (WjV^W^B, 



"K(s) is the perturbed full-order model as in (59) with E 



are as defined in (62) 



(64) 
I, and £ and B 



Corollary 6.2 illustrates that in the case of E = I, both of the reduced system 
matrices, A r and A r , are perturbations of rank min(r, m,p) to the diagonal 
matrix of interpolation points, S. 
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